%% POLICY FUNCTION FOR KNOWN-TYPE CONTRACT
function  [C, Beta, Eta, K, f, ExRet, mu_v] =type_optimal_policy(V0, C0, Beta0, Eta0, params, grid_params)
%% parameters
temp = num2cell(params);
[rho, r, ra, alpha, sigma, lambda, M, eta_min, eta_max] = temp{:};

[n_grid, Y_mat, PHI_mat, dy, dphi] = grid_params{:};

PHI = PHI_mat(1,:)';

%% Value Function
V = V0;

%% Initialize controls
C = zeros(n_grid,1);
Beta = zeros(n_grid,1);
Eta = zeros(n_grid, 1);


options = optimoptions('patternsearch',...
            'UseVectorized',true,'Display','off',...
            'ConstraintTolerance', 1e-12,...
            'FunctionTolerance', 1e-12,...
            'MeshTolerance', 1e-12,...
            'MaxIterations', 1000);

%% Find optimizers
parfor iii = 1:n_grid
    x = PHI(iii);
    v = V(iii);
    
    c0 = C0(iii);
    beta0 = Beta0(iii);
    eta0 = Eta0(iii);

    Aeq = [];
    Beq = [];

    obj = @(z) type_fun(z, x, v, params);

    zz = patternsearch(obj,[c0, beta0, eta0],[],[],Aeq,Beq,[0,0,eta_max],[M,M,eta_max],[],options)

    C(iii) = zz(1);
    Beta(iii) = zz(2);
    Eta(iii) = zz(3);

end

%% Update controls
ExRet = sigma*Eta.*PHI;

K = (Beta.*C.^rho./(sigma*lambda));

mu_v = ra - (C.^(1-rho))/(1-rho) + 0.5*rho*Beta.^2;

f = C - ExRet.*K;
end

